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We present a new fluid-dynamical model of traffic flow. This model generalizes the model of 
Aw and Rascle [SIAM J. Appl. Math. 60 916-938] and Greenberg [SIAM J. Appl. Math 62 729- 
745] by prescribing a more general source term to the velocity equation. This source term can be 
physically motivated by experimental data, when taking into account relaxation and reaction time. 
In particular, the new model has a (linearly) unstable regime as observed in traffic dynamics. We 
develop a numerical code, which solves the corresponding system of balance laws. Applying our code 
to a wide variety of initial data, we find the observed inverse- A shape of the fundamental diagram 
of traffic flow. 

PACS numbers: 89.40.Bb, 05.10.-a, 47.20. Cq 



I. INTRODUCTION 



the ARG model can be rewritten as 



After two-equation models of traffic flow were seriously 
criticized by Daganzo Q the main focus of the traffic 
community has shifted towards microscopic models of 
traffic flow. However, the criticism has been overcome, 
see e.g. 0, Q- By replacing the space derivate in old 
two-equation models by the convective derivative, Aw 
and Rascle and Greenberg |j| deduced a two-equation 
model, which solves all inconsistencies of the earlier mod- 
els as they showed with a detailed mathematical analysis 
and numerical simulations. In particular, in their model 
(in the following called ARG model), no information 
travels faster than the vehicle velocity , i.e. in general 
drivers do not react to the traffic situation behind them. 
Moreover, the velocity is always non-negative. In the 
ARG model traffic flow is prescribed by the following sys- 
tem of balance laws determining the density p = p(t, x) 
and velocity v = v(t, x) of cars 



dp d(pv) 



= 0, 



dt dx 

d(p{v - u(p)) d{pv(v - u(p))) p(u(p) - v) 



dt 



dx 



T 



(1) 
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As usual, (t, x) denote the time and space variable. u(p) 
denotes the equilibrium velocity, which fulfills the follow- 
ing conditions 



u\p) < for < p < p m , 



d 2 ( P u(p)) 
dp 2 



< for < p < p m , 



(3) 
(4) 



with the maximum vehicle density p m . T > is an addi- 
tional parameter, the relaxation time. In the formal limit 
T — > the ARG model reduces to the classic Lighthill- 
Whitham- Richards model la, HOI- For smooth solutions, 



dp dp dv 
dt dx P dx 
dv , . . , , dv 

!Tt +{v + pu{p)) dx 



0, (5) 
u(p) — v 



T 
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In our opinion the ARG model still has a drawback, 
i.e. it can not explain the growth of structures and the 
general behavior for congested traffic, as observed in traf- 
fic dynamics (see e.g. |9l Il0|). To see this, we consider a 
linear stability analysis around the equilibrium solution 
p(t, x) = p Q , v(t, x) = u(po), i.e. 

p(t,x) — po + penp(ikx + u>(k)t), (7) 
v(t,x) — u(pq) + v exp(ikx + ui(k)t). (8) 

Substituting this ansatz into system Q-lj^J we obtain 



iku ikpo 
f uj + ± + ik(u + p u') J \v 



(9) 



Nontrivial solutions of this linear system exist if and only 
if 



(u + iku) (uj + — + ik(u + pou'ij + i 

or equivalently for 

uii = -ik(u + p u r ), 
ix>2 = — ^ — iku. 



*^ = 0, (10) 



(11) 
(12) 



For the stability properties, the real parts of the above 
solutions are important, i.e. 



i?e(w 2 ) 



0. 



1 

T' 



(13) 
(14) 



For T > both real parts are nonpositive, which means 
that the ARG model is linearly stable, the velocity v re- 
laxes to the equilibrium velocity u in the entire region 



2 



< p < p m - This is clearly in contrast to observations, 
where a wide range of states in the fundamental dia- 
gram, the relation between vehicle flux and the density, 
are observed for congested traffic flow. To cure this de- 
fect, Greenberg, Klar and Rascle developed an extended 
model with two equilibrium velocities In the paper 
here, we propose an alternative model, which takes into 
account the reaction time of drivers (as well as mechani- 
cal restrictions). 

We give a physical argument for our new model and 
define it in Sec. [HI Section lTTTI presents the methods used 
for numerically solving the model equations. Section llVI 
describes tests to validate our numerical algorithm, be- 
fore we finally discuss the numerical results on the fun- 
damental diagram obtained with our model in Sec. 



II. A HEURISTIC DERIVATION OF THE NEW 
MODEL 

Before we turn to the new model, let us first give a sim- 
ple derivation of the ARG model. Note, that the model 
was mathematically derived from car following theory 
in |l2j | . Suppose that in the reference frame of individual 
drivers, drivers adjust their speed v in such a way, that 
they asymptotically approach the equilibrium velocity u, 
i.e. 



d(v — u) u 



dt 



T 



(15) 



Here, T = const > is the relaxation time. In compar- 
ison to optimal velocity models (see e.g. 0) the equi- 
librium velocity term on the left has been added which 
vanishes for u = const, ft is easy to verify, that the ana- 
lytical solution of the ordinary differential equation i|15l) 
reads 



u(t)=«(t) + (i;(0)-«(0))ex P (--) 



(16) 



In the coordinate system of the road, Eq. I|15f) translates 
to 



dt 

Moreover, since 

du\ 



d(v — u) | d[v — u) u — v 



dx 



T 
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pu 



, dv 
dx ' 



(17) 



(18) 



where we have used the continuity equation JSJ for the 
last equality, we recover the velocity equation of the ARG 
model © ■ From this derivation, it is obvious that drivers 
instantaneously react to the current traffic situation. 

We therefore tried to generalize Eq. (|15fl and took the 
reaction time of drivers r > into account 



dv 
~dt 



(M) 



du 
~dl 

u(p(x - 



(p(x — VT, t 
- VT, t — t)) - 



•r)) = 
v(x 



(19) 



VT, t 



Using a Taylor series expansion in r and keeping only 
terms up to order in r and T, i.e. 



du du(p{x,t)) du{p{x,t)) 

- {p ( X -VT,t-T)) = - +V ^ 



+0\t,T), 

u(p(x — VT,t — t)) = u(p{x,t)) 



(20) 



-™'(p(x, 
+0 2 {t, T) 

= u + pu'^T + 2 (T,T), (21) 
dx 

i \ , x ( dv dv 

v{x~VT,t-T) = V(x,t) - T — + V — 

\ dt dx 



we find 



+0 2 (t, T), 



dv , .. dv u(p) 



(22) 



(23) 



This equation is identical to the velocity equation of the 
ARG model ©, except that the relaxation time T has 
been replaced by T — r. In particular it follows from the 
stability analysis of the ARG model, that for r > T the 
new system is (linearly) unstable. 

Before we look at the experimental data on the re- 
laxation and reaction time, we remark that one could be 
tempted to include an anticipation length into the model, 
as e.g. in [T3. ITfij. This approach has not been followed 
here for two reasons: First, the ARG model already in- 
cludes anticipatory elements, as noted by 0- Second, 
including the anticipation length into the above deriva- 
tion yields a system, which does not guarantee that the 
maximum speed at which information travels is bounded 
from above by the velocity of cars, and is therefore unre- 
alistic. 

For the reaction time r, typical values are of the order 



0.5 .. 1 s. 



(24) 



Fig-fflshows experimental results for the relaxation time 
f taken from the review article |l6j| . However, these val- 
ues have to be interpreted with care and cannot be trans- 
lated directly to our model context. To see this, we note 
that the relaxation time T is determined for the ansatz 
v(t,x) = u(p(t+T, x+Ax)), i.e. after the relaxation time 
T the driver has fully adjusted to the equilibrium velocity 
u. Here, according to Eqs. (|15fl and (|16fl . the equilibrium 
velocity in general will never be reached exactly. Instead, 
if we require, that \v(t) - u(t)\ < \v(0) - u(0)|/1000, we 
find that t > 6.908 T rs T. I hence seems reasonable to 
set 



T : 



10 T. 



(25) 



Taking average values r = 0.75 s and t — 7.5 T, we 
indeed find that T — r < for about p > 50 [1/km/lane]. 
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FIG. 1: Dependence of the relaxation time T on the vehicle 
density per lane. Taken from ITr^| . 



It was also pointed out in ^(| that for large densities the 
relaxation time T increases, which means that T — r > 
for p k p m . 

One could try to repeat the derivation leading to 
Eq. (|23|) for a general relaxation time T = T(p, v). Note, 
that the above derivation is only valid for a constant 
relaxation time. Moreover, it involves only the leading 
term of a Taylor series expansion. We therefore decided 
to generalize the velocity equation of the ARG model in 
the following way 

^ + (v + pu'(p))^ = (3(p, v)(u(p) - v). (26) 

Note that we do not require v < u as Greenberg 
From the experimental data and the argument put for- 
ward before (note that the sign of f3 determines whether 
the traffic flow is linearly stable or not) we require 

/3(p, v) < for < p x < p < p 2 < p m ,v = u(p), (27) 

(28) 



lim /3(p,v) > 0, 
lim f3(p, v) > 0. 



(29) 



and the parameters p m = 300 [1/km], u m = 140 km/h, 
ri\ = 0.35, 7i 2 = 1 (note that with these parameters, 
the equilibrium velocity of Cremer 1)32(1 fulfills the condi- 
tions © and gj), f = 1 s, a = 12, p x = 70 [1/km] and 
P2 = 270 [1/km] the function P(p,v) already fulfills all 
requirements l(27|l - l(29(l . However, due to mechanical re- 
strictions, the maximum acceleration a c and deceleration 
d c give stronger limitations, i.e. 



dv dv 

—— < a c and — > d c 

dt ~ dt ~ 



(33) 



with typical values a c = 2 m/s 2 and d c = —5 m/s 2 . Since 
the resulting system is not strictly hyperbolic for equality 
in Eq. I(33|l . which is problematic for a numerical solution, 
we prescribe the limitations on 



d(v — u) , d(v — u) 

' < a r and ■ > d r 



dt 



dt 



(34) 



which then leads to the functional form (|30l) . We plot the 
function (3{p, v) for the mentioned parameter values in 
Fig. [3 We stress that the above functions describe reality 
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Throughout this paper we use a functional form 

^, X0(p,v)(u-v)-a e >O, 

P(P> v ) = { #5. if Hp, v)(u -v)-dc<0, 
j3(p,v), else, 



where the function /3(p, v) is defined as 
TV 



(30) 



= " 1 1 + J^A + J-(-( Pl + p 2 )p + p 2 ) ) . ( :-! ! , 
T V u m pip 2 



For the choice of the velocity-density relation of Cre- 
mer [13 



(32) 



FIG. 2: Function /3(p,v) defined in Eq. l(3u). We used units 
[p] = 1/km and [v]=km/h and [f3] = 1/h. 

only qualitatively. For realistic simulations of traffic flow, 
experimental data are required to determine (3{p, v). 

In the new model, traffic flow is described by the fol- 
lowing system of balance laws 



dp d(pv) 
dt dx 
d(p(v~u(p)) d{pv{v - u{p))) 



dt 



dx 



0, (35) 



or equivalently for smooth solutions by Eqs. (0) and J2EJ- 
As the corresponding system of the ARG model, the new 
system is strictly hyperbolic for < p < p m . 
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III. THE NUMERICAL IMPLEMENTATION 



Explicitly, 



Writing traffic flow as a system of balance laws in 
Eqs. i|35|l and i|36|) is very adequate for numerical pur- 
poses, as it allows the application of well-established hy- 
drodynamic methods for the numerical solution. We use 
a high-resolution shock-capturing scheme with approx- 
imate Riemann solver for the numerical solution (see 
e.g. 0). 

We rewrite Eqs. (|35|l and (|36(l in the form 



dU 
~dt 



dF(U) 
dx 



= S(U), 



where 



U = 



F{U) 



S{U) 



P 

p(v — u) 
pv 

pv(v — u] 


f3p(u — v 




t/i«(£/i) N 
U 2 u(U 1 ) i 



(37) 

(38) 
(39) 
(40) 



We use the second-order reconstruction scheme of van 
Leer to reconstruct quantities at cell interfaces. At 
cell i with cell center at the location = xq + iAx the 
update in time from t n to t n+1 is performed according to 
an conservative algorithm 



UT +1 = m 



,) + AtS l 



(41) 



- t Tl 



To obtain 



where [/" = U(x t ,t n ) and At = t r ' 
a higher order of convergence, we use the third order 
scheme of Shu and Osher [2(j]. The numerical fluxes 
F are determined according to the flux-formula of Mar- 
quina 21], which reads 



F 



1 



F 1 



Aq). 



(42) 



Here, the superscripts R and L denote the reconstructed 
values on the right and left of a cell interface. The nu- 
merical viscosity term takes the form 



Aq = R^IAUaxLV - R L \A\ max L L U L . (43) 
The matrix |A| max involves the characteristic speeds, 



= ^maxdAfUAfl! 







max(|A»|,|A^|) 

where the characteristic speeds read explicitly 

Ai = v + pu', 
A 2 = v. 



(44) 



(45) 
(46) 



R and L are the matrices of the right and left eigenvectors 
of the matrix 



dF 
dU 



pu 



1 



— (v — u) 2 + p(v — u)u' 2v — u 



(47) 



R 



u ~ v v — u — pu 



k 1 f u — v + pu' 1 



pu 



v — u 



(48) 
(49) 



IV. CODE TESTS 

We checked, that our numerical algorithm is conver- 
gent. Moreover, the density equation l|35|l is a strict con- 
servation law. Prescribing periodic boundary conditions 
as in Sec. ]V\ the total number of cars included in the 
numerical domain il should therefore be constant, i.e. 



p dx = const. 



(50) 



We checked, that our numerical code fulfills Eq. (|5*0|l up 
to machine precision (see also the corresponding results 
for a network simulation based on the Lighthill-Whitham 
theory in [23]). 
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FIG. 3: Numerical solution for a Riemann problem of Aw and 
Rascle 0- The numerical solution at time t=0.03 h (solid 
line) for the initial data (dashed line) reproduces the exact 
solution (dotted line). 

Finally, Aw and Rascle presented in their paper the 
exact solution of a Riemann problem, for which old two- 
equation models fail to describe the correct behavior 
(see Q, Fig. 5.4). This Riemann problem consists of 
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FIG. 4: Time evolution of the density for stable and insta- 
ble initial data. We prescribe an equilibrium density po = 65 
[1/km] and po — 80 [1/km] respectively, and on top a sinu- 
soidal density perturbation. For the instable data, the initial 
perturbation located at x = 2.5 [km] is amplified and finally 
two clusters form. 

the following initial data 

_ J 0, if x < 1 km, , . 

9 ~ \ p+, if x > 1 km, 

v = (°' ^tlt™' (52) 
I it x > 1 km. v ' 

The exact solution to the homogeneous system consists 
of the constant state on the right (p + , v+) moving to the 
right with velocity v + , leaving behind vacuum. If we 
chose v + — u(p+), this exact solution will carry over to 
our inhomogeneous system. Fig. [21 displays our numer- 
ical solution for a choice p+ = 50 [1/km]. For numeri- 
cal reasons, we prescribe a density p = 10~ 6 [1/km] for 
x < 1 km. Note, that our numerical algorithm resolves 
the steep gradient within only a few grid cells, at the 
same time reproducing the correct velocity at which the 
constant state moves to the right. Moreover, the velocity 
relaxes to the equilibrium velocity behind the constant 
state (p+, v+). 

V. RESULTS ON THE FUNDAMENTAL 
DIAGRAM 

For the results presented in this section we restrict us 
to a 7 km long section of a (two-lane) highway with peri- 
odic boundary conditions. On this section of the highway, 
we start our simulations with constant equilibrium data 




FIG. 5: Time evolution of the velocity for stable (po = 65 
[1/km]) and instable initial data (po = 80 [1/km]) with initial 
density perturbation. 

p = pa, v — u(po) and in addition between kilometers 2 
and 3 a sinusoidal density perturbation 

Ap = sin(7rx) for 2 < x < 3 km. (53) 

For all presented numerical results we used a resolution 
of 50 m. Fig. 01 shows the evolution of these data for 
parameters po — 65 [1 /km] and po — 80 [1 /km] . Whereas 
the amplitude of the perturbation is gradually damped 
with time for the stable initial data po = 65 [1/km], the 
amplitude of the perturbation increases for the instable 
initial data po = 80 [1/km]. Moreover, the perturbation 
travels with a larger velocity downstream in the first case. 
For the unstable situation, two clusters are forming. We 
plot the corresponding time evolutions of the velocity in 
Fig. 

To obtain a more general picture we varied the initial 
density in the entire density regime and analyzed the re- 
sulting flow-density-relation as a function of time. More 
precisely, we used initial values for the equilibrium data 
po = 2,4, ....,298 and read off the resulting values for 
the density p and the flux function pv at 5 equidistantly 
distributed cross sections of the highway. Fig. shows 
the results for evolution times t = 0.00 h (initial data), 
t = 0.05 h, t = 0.10 h, t = 0.15 h, t = 0.20 h and 
t = 0.25 h. For the initial data, the flow-density curve 
corresponds closely to the equilibrium flow density, the 
initial perturbation (|53|l being negligible for the visual 
output. After an evolution time t = 0.05 h, the equilib- 
rium flow-density-curve is still visible, but in the unstable 
regime for densities 70 < p < 270 [1 /km] two new flow- 
density-curves start to appear. In the evolution further 
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the two new branches produce an inverse- A shape. 



VI. CONCLUSION AND OUTLOOK 

We generalized the traffic model of Aw, Rascle and 
Greenberg by prescribing a more general source term 
to the velocity equation and developed a new numeri- 
cal code to solve the resulting system of balance laws. In 
total our (numerical) results show: 
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• The new model can explain the large variance of 
the measured values of the fundamental diagram in 
the congested regime, which correspond to fluctua- 
tions between two branches in the unstable density 
regime. Moreover, due to the stability properties, 
the model predicts oscillations in the relative ve- 
locity of cars in the congested regime, as they are 
found in experimental data. At the same time, it 
reproduces the small variance of velocities for free 
traffic flow and can explain the appearance of wide 
traffic jams. 
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FIG. 6: Fundamental diagram for the initial data (t = 0.00 h) 
(first row on the left), for t = 0.05 h (first row on the right), for 
t = 0.10 h (second row on the left), for t = 0.15 h (second row 
on the right), for t = 0.20 h (third row on the left) and for t = 
0.25 h (third row on the right). At intermediate densities, our 
traffic model is linearly unstable, the representative points in 
the fundamental diagram are shifted towards two branches, 
which gives the visual impression of an inverted A. 



in time, the equilibrium density curve vanishes. Instead 



• Macroscopic traffic models have often used an equi- 
librium velocity u(p), for which > in the 
congested regime, in order to account for the val- 
ues of traffic flow at the maximum (the tip of the 
inverted A). According to our study, this is not 
necessary, as the high values for the fluxes can be 
explained with overcritical solutions and an equi- 



librium velocity function with 
where. 



d(pu(p)) 



< every- 



The new model, which is a deterministic and effective 
one-lane model, has the capacity of reproducing many 
features observed in traffic dynamics. In the presented 
work, the form of the function (3 in Fig. |3 was motivated 
by a physical argument, but the quantitative details were 
determined rather ad hoc. However, we found that the 
fundamental diagram in the unstable region (e.g. the tip 
of the inverted A) depends on the particular form of /3. 
Hence one should try to determine the function (3 from 
experimental data of the fundamental diagram. In our 
opinion, the presented algorithm is adequate for the use 
in network simulations of traffic flow. 
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